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Synchronization of coupled oscillators on a d-dimensional lattice with the power-law coupling 
G{r) = go/r" and randomly distributed intrinsic frequency is analyzed. A systematic perturbation 
theory is developed to calculate the order parameter profile and correlation functions in powers of 
e = a/d — 1. For ct < d, the system exhibits a sharp synchronization transition as described by the 
conventional mean-field theory. For a > d, the transition is smeared by the quenched disorder, and 
the macroscopic order parameter -ip decays slowly with go as \ip\ oc (?o- 
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Introduction. Collective oscillations of active inter- 
acting elements are observed in a variety of physical, 
chemical, and biological systems far from equilibrium. 
Numerous studies have been devoted to the mutual en- 
trainment of oscillators that have different intrinsic fre- 
quencies [H-ll]. A class of models with global (or mean- 
field) coupling have enjoyed deep theoretical understand- 
ing y, Q . The phase of the oscillators become coherent 
as the coupling strength go exceeds a threshold, which 
is the onset of synchronization. Extensive research has 
been focused on the transition behavior The am- 

plitude of the macroscopic order parameter scales as 
IV"! oc (50 — gc)^ , with (3 = 1/2 for the original mean- 
field model by Kuramoto [1], while /3 = 1 for some other 
types of coupling 

Compared to the case of global coupling, behaviors 
of locally [1, [13] and non-locally 11, [11] coupled os- 
cillators are still widely open questions. In particular, 
knowledge about synchronization caused by long-range 
interactions is quite limited (l3l - [l5| . although they are 
ubiquitous in Nature in the form of, e.g., gravitational, 
electromagnetic, elastic, and hydrodynamic forces. Early 
numerical works for the power-law coupling cx 1/r" in 
d-dimensional array of oscillators show that glob al syn- 
chronization is possible for a < 2 {d ~ I) [13|, while 
system-size effect is significant for a < d {d — 1, 2) [l3 |. 

Recently, we proposed a simple model of microfluidic 
carpets [ig, [l3| , which is a two-dimensional array of ro- 
tors with a hydrodynamic coupling oc 1/r^. The model 
exhibits an unconventionally smooth transition to the 
synchronized state [17j . The macroscopic order param- 
eter decays gradually as the randomness is increased, in 
contrast to the sharp transition for global coupling. 

Motivated by the numerical results, this Letter theo- 
retically addresses synchronization of oscillators with a 
general class of long-range coupling. We will develop 
a systematic perturbation expansion around the mean 
field, taking the moments <t„ of the interaction G{r) as 
the small parameters (which is analogous in spirit to the 
cluster expansion in the classical gas theory). For the 
power-law coupling G(r) = g^/r"', it is equivalent to a se- 
ries expansion ine^a/d— 1. We will solve for the order 



parameter profile and correlation functions up to 0{e^). 
The main finding of this paper will be that the macro- 
scopic order parameter for a > d behaves as oi g^ for 
5o — > 0, which means that synchronization persists for 
arbitrary weak coupling. We interpret it as the result of 
quenched spatial heterogeneity. In contrast, for a < d, 
the heterogeneity is averaged out and the transition is 
exactly described by the mean-field theory. 

Model. In our model, oscillators indexed hy i = 
1,2, . . . , N are arrayed on a d-dimensional regular lat- 
tice with the unit grid size. The phase 0i of the i-th 
oscillator located at obeys the dynamic equation. 



dt 



= u},-Y^ G{r, - Tj) sin {<j), - , (1) 

where uJi is the intrinsic frequency that has the Gaussian 
distribution with the standard deviation ujq, 
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We require the coupling function G(r) to be positive, 
slowly decreasing function of ]r], so that its moments 



G(r, 



rapidly decays with n. To be specific, let us consider 
the power-law coupling G(r) = gojr'^ with the constants 
go > and a > 0. We normalize the coupling by rescaling 
time so that ci = 1 without losing generality. For the 
global coupling (a = 0), we have G(r) = g^ — l/N, and 
the moments cr„ = 1/A^"^^ for n > 2 vanish as — ?> 00. 
In more general, for a < d, the integral / d'^r/r°' diverges 
with the system dimension Tmax ~ N^^'^, which means 
that go ~ N"/'^^^ and cr„ (n > 2) vanish as N ^ 00. 
This is true also for a = d, except that the divergence of 
go is logarithmic. On the other hand, for a > d, we have 
go < 1 and tT„ « jiq — >■ (n — )■ 00). Regarding e = a/d—1 
as the small parameter, we can show that (t„ — 0(e"). 
For example, for d = 1, we have go = 1/2^(1 -I- e) ~ e/27, 
and (T„ = 2({na)go ~ 2({n){e/2j)'"' for n > 2, where 7 is 
Euler's constant. Our perturbation theory will be given 
as a series expansion in e via cr„ = 0(e"). 
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Order Parameter. In order to describe the collective 
behavior, we introduce the site-dependent complex order 
parameter 'ipi with its amplitude pi and phase 9i 



with which we can rewrite ([TJ as 



dt 



- Pi sm(0i 



(2) 



(3) 



Note that pi < 1 due to the normalization of ci. When 
the coupling is long-ranged, "0^ involves infinitely many 
oscillators and is expected to change much slower than 
(pi. Therefore, we approximate V'i to be constant in time. 
Then Eq.([5]) is replaced by its temporal average. 



(4) 



where — G{ri — rj) for i j, and Gij = for 
i = j. The function E(pj,ujj) is the temporal average 
of e^^'^^~^^\ and is calculated following the original pre- 
scription by Kuramoto 0, 13]- First, for an oscillator that 
satisfies \uji\ < pi (coherent case), Eq.Q allows the sta- 
tionary solution (pi = 9i + sin~^ (uji/pi), which gives 



^1 -^i 
-k- + « — . 



On the other hand, if \LUi\ > pi (incoherent case), 
Eq.(|31) has a drifting solution, which visits each value 
of (pi with the frequency that is inversely propor- 
tional to the angular velocity: v{(j)i) = Vi\(j)i\^^ = 
Vi\uJi — PiSm{(pi — 9i)\ . Here, the constant Vi — 
■^yj^'i ~ Pi ensures the normalization J^^ dpi'{p) — 1. 
It gives the temporal average J^^ d(j}iy{pi)e^^'^^^^^^ as 



E{pj,ujj) 



. Wo 



Pj 



Perturbation Expansion. Let us introduce the two- 
dimensional vector t/jj = {'4'iRj'4'ii) — {pi cos 9i, Pi sin6'i), 
to rewrite Eq.(|4|) in the vectorial form 



F{rp„^j) 



Re 
Im 



,^'^E{p,,u:,), 



(5) 



(6) 



Our task is to calculate the spatial average of the order 
parameter, 

^a^J^'^i^^a, a = R,I, 



which is equivalent to the ensemble average over w^'s. 
Expecting that spatial fluctuation of the order parameter 
is small for long-range interactions, we expand the RHS 
of ([5]) with respect to the deviation 5%pj = xp ^ — ip, as 



Fja.bc 



SipjbSipjc 



,(7) 



with Fj, 

F. 



= Fa{ip,ujj), Fja,b = ^^Fjai'ip,ujj), and 

-Fja{ip,ujj). Here and hereafter, summa- 
R,I and i,j, k,£ = 



tion over repeated indices a, b, a, d 
1,2,...,N are implied. We decompose the zeroth and 
first order coefficients into their averages fa = /a('0) = 



(Fja)^. , 
^Fja.b - 



fa,b 



and the deviations 5Ka 



/a, I 



Subtracting ■0^ from ([7]) and 



then multiplying by the inverse of the 2N x 2N matrix 
^Kib = ^ij^ab - Gijfa,b, which is expanded as 
SijSab + Gijfa,b + Gijf^fj + ... with Gij — GikGkj, 
fab = fa.cfc.b, etc., we obtain 



5V„ = Aa+r^,Ujb, 

Ujb = SFjb + 5Fjb,cSijjjc 

df 



-Fjb,cd6i'jc5'Pjd - 



ah 



(8) 
(9) 

(10) 



r^, = [M-\G®h)X^^ = G,,5ab + Glfa^b + ...,m 

where we used the 2x2 matrices I2 — {5ab} f^nd 
1= = {fa.b}- Eas. (|8l9l) can be diagrammatized as shown 
in Figlljb), by combining the symbols defined in Figllja). 
Recursively using ^ for the Sxp's in ©, we get an ex- 
pansion of 5i)ia in terms of A^, T^^^, SFjaitp), FjaMi^P) 
and their derivatives; see FiglUc). The terminators Aq 
and SFja are connected by the vertices SFja,b, Fja.bc, ■ ■ ■ 
and propagators F^'^^ to the site i. For example, the graph 
framed by solid lines reads F^-^j^ • SFjb^c ■ ^^j^SEkd, and the 
dot-framed graph reads T'^^ ■ \F,b,cd ■ ^SFke ■ T'f^SFif. 

Now we take the average of Eq.® over the distri- 
bution of Wi's. The LHS vanishes by the definition of 
dip. On the RHS, SFjb and its derivative SFjb_c are av- 
eraged out unless they are correlated with a partner at 
the same site. Graphically, it means that the legs of the 
graphs (with the black dots at their ends) have to be at- 
tached to each other or to vertices to produce correlation 
terms. For example, the dot-framed graph in Fig[ljc) 
yields the corresponding graph in FigUJd), which reads 
^'ab ■ 1 (Fjb.cd)^ ■ n^eK^ (SFkJFkf),, where (• • means 
the average over the distribution P{ujj). Using the ex- 
pansion PT|) with the trace G"^ — an, we obtain the 
O(e^) expression of this graph as 



0-2 



fa,cd {dcd fcfd) 1 



where the functions fa,cd{ip) = gU^^ and gab{ip) = 
{FkaFkb)k ^'''6 introduced. Another graph that gives an 
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X A,, 



(a) ° 5v|f„ 

(b) □ = X 

SVia A 3 Tab 5Fjb r,b 5Fjb ,.8v|/jc 



-ex + — ® + 



(C) □ = X + • + ■ 

+ — o O x + — e— ^ + 1 — sT"* I 

X — • . 



■e- 



(d) = X + 



FIG. 1: Diagrammatic expression of the perturbation scheme, 
(a) Definition of symbols, (b) Expression of Eas. (|8l9l) and (c) 
its recursive expansion, (d) The ensemble average of (c). See 
text for interpretation of the framed graphs. 



O(e^) contribution is framed by solid lines in Fig[ljd). It 

reads ^ab^cd {^Pjb,cSFjd) ^ and is approximated as 

0'2fc,d {hda,c " fdfa,c) 

with the function hab,c{tl>) = {FjaFjb,c)j- We can see 
that these two graphs and are the only O(e^) contri- 
butions. Combining them and using Eg. pH)) . we obtain 
the self-consistent equation for if) to O(e^) as 



i'a = fa + Cr2{Sab ~ fa,b)Vb, 



(12) 



1 



Vb - fcAhdb,c-hh,c) + -^fb.cd{gcd-fch).{li) 

Correlation Function. The correlation function of the 
order parameter C*^{, = Caf,(ri — r^) = (Stl^iaS^jb) can be 
also computed using the diagrams. There is only one 
non- vanishing graph at O(e^), which gives 



'^ab - ^ij i9ab - fafb) ■ 



(14) 



Note that Gf^ — GikGkj is a function of r.y = — r^. At 
large distance, it decays as Gf^ cx |rij|~''(^+^'^) for e > 0, 
as we can see from a simple dimensional analysis. (For 
d>2, Gfj depends also on the direction of r^j reflecting 
the lattice anisotropy). On the other hand, setting i — j 
in (jl4p . we obtain the variance of the order parameter. 



(SlptaS'lptb) = 0-2 (gab - fafb) 



(15) 



Transition Behavior. In order to solve the self- 
consistent equation p2ll3p . we need to compute fa, gab, 
hab,c and their derivatives as functions of V' = pe'^. To 
simplify calculations, we choose the coordinate frame in 
which = 0. Then the ensemble average of Eq.® gives 



fa = ea(p) 



dujP{uj)Eaip,u}) 



(16) 



(a = R,I). Here, Er{p,uj) and Ej{p,uj) are the real 
and imaginary parts of E{p,uj), respectively. Note that 
// = e/(p) = thanks to the parity of P{oj) (even) and 
Ej{p,io) (odd). The quadratic moments read 

9RR = eRRip), gRi = giR = 0, gii = eii{p), 

eabip) = / dujP{uj)Ea{p,uj)Ei,(p,uj). (17) 



The calculations of the derivatives fa^b, fa.bc and hab,c are 
also straightforward. The non-vanishing components are 



- e'^, //,/ — sr, fn^RR 
e^, and Hrr^r = e'j^j^ 



,/2, h 



RIJ 



found to be fR^R 

fl.RI = fl.IR = 

eRR, hiRj = -e//, hii,R = e'jj/2, where ' = d/dp and 
the abbreviations br — cr/'p, crr, = eRRj'p, and en = 
eji/'p are used. Substituting these into Eqs. p2ll3p . we 
obtain 



'4'R 

Vr 



sr 

1 

2 



a2{l-e'^)VR, 
e-R {e-RR - 2ei?e^) - 2?fl,e// 
+ eR {eRR - e\) + ^^e// 



(18) 



(19) 



and if) J — 0. On the RHS of (fT8)) are functions of p, which 
is related to -0^ = pcosO on the LHS via the expansion 

^R + {5i>'j)/l^R + O(e^). Using this 



in the RHS of ([T5|) with the result (Stpj) — (T2eii taken 
from Eg. (fist , we arrive at the final form of the O(e^) 
self-consistent equation. 



V'fl = efl, (72 [(1 - e^)VR + e^e//] , 



(20) 



with the functions on the RHS evaluated at p = ipR. Its 
solution gives the order parameter profile ^r = 'iPr{ujo). 
For f72 = 0, or a < rf, Eq. d^Dl) reduces to the mean- 
field equation = eR{ipR) tL,!^]. The Taylor expansion 

BRip) ~ {uJc/uJo){p - p^/8) with UJc = 



0.627 

reproduces the global-coupling result that the order pa- 
rameter vanishes for ujq > Uc- In contrast, for (72 > 0, or 
a > d, there is no sharp transition, and the order param- 
eter exhibits a long tail at large wq, In fact, the approx- 
imation eRR{p) ~ eii{p) ~ {8u!c/3TruJo)p for p ^ 1 gives 
the asymptotic behavior of -0^ for wo ^ Wc, 



(72 

6w2- 



(21) 



The complete order parameter profile is obtained by nu- 
merical computation of the functions ea(p) and eab(p), 
and is shown in Figl^^a). Note that ipR in the current 
coordinate frame corresponds to \ip\ in the general frame. 
As we can see, the deviation from the mean-field profile 
is significant even for relatively small values of CT2 . (For 
comparison, (T2 = 0.2 for {d,a) — (1,2) and (72 — 0.057 
for {d,a) = (2,3) (square lattice).) The macroscopic 
order parameter is larger than the mean-field value for 
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Wo > Wt ~ 0.504, and smaller for wq < oJt for any non- 
zero value of (72. The enhancement of synchronization 
for large Wq might look counter-intuitive, but it is a nat- 
ural result of the spatial heterogeneity; there are regions 
that are more uniform than the others in terms of the in- 
trinsic frequencies of the oscillators they contain. These 
regions can remain synchronized when the other regions 
are desynchronized, and contribute to the long tail of the 
order parameter profile. This effect of quenched hetero- 
geneity is averaged out in the global-coupling case. Note 
also that we have rescaled the timescale so that cti = 1. 
If <Ji is not normalized, we must divide the intrinsic fre- 
quency and the order parameter by (Ti, which modifies 
Eq. ([2T|) as IV'I ~ cricr2/6aJo = C'ffo/'^O' where C is a func- 
tion of a and d. 



(a) 



1^ 




0.2 0.4 0.6 0.! 



(b) 



> 

CO 



0.5 
0.4 
0.3 
0.2 
0.1 




02=0-00 — 

0.01 ■ ■ 

^ 0.03 - - 

"N 0.10 -- 
/ \ 0.30 



0.2 0.4 0.6 0.! 



1 1.2 1.4 



FIG. 2: (a) The macroscopic order parameter and (b) 
the standard deviation std{ip) = (\5^\ 



2\l/2 



as functions of 



ujo. The long tails scale as oc a^jio^ and std(?/;) oc a2/ojo. 

It should be briefly mentioned that the self-consistent 
solution bifurcates at very small into two stable 
branches — 1 ^^^^ ?/'fl2 ^ 1- The threshold uJb rises 
with (72; e.g., ujt — 0.003 for a2 = 0.1 and uJb = 0.022 
for (72 = 0.3. However, it turns out that the lower 
branch docs not satisfy the condition for the series (fTTj) 
to converge. It converges with its trace X)^o '''"/"a if 
max(e^(V'i^), ei?(i/'^)) < l/go- Plotted in FiglSlJa) is the 
upper branch, which always meets the condition. 



1/2 

The standard deviation std(?/') — (jf^i/'p) is readily 
calculated from Eqs. (|15ll6ll7p . and is plotted in FiglSJ^b). 
For any non-zero value of cr2, it exhibits a peak near 
LdQ — ujc and a long-tail for ujq 3> Wc- The asymptotic 
behavior for large ujq is obtained via the Taylor expansion 
of e_R(p), eflfl(p) and e//(p), as std{i>) w a2/^^/oJ^ujo- 



Summary. We have found that the mean-field pic- 
ture of sharp synchronization transition is valid only for 
a < d, and the transition is broadened for a > d. It 
could be regarded as a novel example of smeared transi- 
tion in random systems, which usually requires spatially 
correlated disorder Is'] . The limitations of the perturba- 
tion theory for large a should be assessed by analysis of 
higher order corrections and comparison with numerical 
results, which are beyond the scope of the present paper 
and will be discussed elsewhere. 

I wish to thank Ramin Golestanian for useful com- 
ments, discussions, and collaborated works that moti- 
vated the present study. 
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